/**
 * @file basic_rng.h
 *
 * Random number generators for basic distributions
 *
 * @author Dahua Lin
 */

#ifdef _MSC_VER
#pragma once
#endif

#ifndef GMLPP_BASIC_RNG_H_
#define GMLPP_BASIC_RNG_H_

#include <gmlpp/base/basic_defs.h>
#include <gmlpp/base/basic_math.h>

namespace gml
{

	template<class RStream>
	class RandomGen
	{
	public:
		typedef RStream rand_stream_type;

	public:
		RandomGen(RStream& rs)
		: m_rs(rs)
		, m_randn_cached(false), m_randn_cv(0)
		{
		}

		GML_ENSURE_INLINE uint32_t randi(uint32_t m) // --> discrete uniform [0, m-1)
		{
			return m_rs.gen_u32() % m;
		}

		GML_ENSURE_INLINE bool randb(double p) // --> bernoulli(p)
		{
			return m_rs.gen_f64() < p;
		}

		GML_ENSURE_INLINE double randu() // --> real uniform (0, 1)
		{
			return m_rs.gen_f64();
		}

		double randn()  // --> normal
		{
			if (m_randn_cached)
			{
				m_randn_cached = false;
				return m_randn_cv;
			}
			else
			{
				double rv;
				marsaglia_polar_randn(rv, m_randn_cv);

				m_randn_cached = true;
				return rv;
			}
		}

		GML_ENSURE_INLINE double rande()  // --> exponential
		{
			return -math::log(m_rs.gen_f64());
		}

		GML_ENSURE_INLINE double randg(double alpha) // --> gamma
		{
			if (alpha >= 1)
			{
				return marsaglia_randg(alpha);
			}
			else
			{
				double y = marsaglia_randg(1.0 + alpha);
				return y * math::pow(randu(), 1.0 / alpha);
			}
		}

	private:
		void marsaglia_polar_randn(double& x, double& y)
		{
			double u, v, s;
			do
			{
				u = 1.0 - 2 * m_rs.gen_f64();
				v = 1.0 - 2 * m_rs.gen_f64();
				s = u * u + v * v;
			}while(s >= 1);

			double r = math::sqrt((-2) * math::log(s) / s);
			x = u * r;
			y = v * r;
		}

		double marsaglia_randg(double a);

	private:
		rand_stream_type m_rs;

		bool m_randn_cached;
		double m_randn_cv;

	}; // end class RandomGen


	template<class RStream>
	double RandomGen<RStream>::marsaglia_randg(double a)
	{
		double d = a - 1.0/3;
		double c = 1.0 / sqrt(9.0 * d);
		double v = 0, x = 0;

		for(;;)
		{
			do
			{
				x = randn();
				v = 1.0 + c * x;
			}while(v <= 0);

			v = v * v * v;
			double u = randu();

			double x2 = x * x;
			if(u < 1.0 - 0.0331 * (x2*x2))
				break;

			if (math::log(u) < 0.5 * x2 + d * (1.0 - v + math::log(v)))
				break;
		}

		return d * v;
	}


}

#endif /* BASIC_RNG_H_ */
